#ifndef _SM_SYSTEM_H_
#define _SM_SYSTEM_H_
#include<string>
#include"help.h"
#include"base.h"
#include"orbital.h"
#include"twobodystate.h"
#include"system_ho.h"
#include"system_hf.h"
#include"system_gammow.h"

using std::string;

namespace SM
{
  typedef complexd DataType;
  // typedef double DataType;

  ///orbitals with l,jj,mm,tz
  class Orbital_SM_Tz_M:public Orbital_ljjmmtz
  {
  public:
    Orbital_SM_Tz_M(){}
  Orbital_SM_Tz_M(int _l,int _jj,int _mm,int _tz,DataType _e):Orbital_ljjmmtz(_l,_jj,_mm,_tz),e(_e){}

    DataType e;
  };

  
  //****************************************************************  
  ///system hold interaction in M scheme
  class System_Tz_SM_M:public System<Orbital_SM_Tz_M,TwoBodyState_TzParM,DataType>
    {
    public:
      //typedef System_Tz_Gammow VSystem;
      typedef System_Tz_HF VSystem;
      //typedef System_Tz_HO VSystem;//interaction system
      typedef Orbital_SM_Tz_M OrbitalType;
      typedef TwoBodyState_TzParM TwoBodyStateType;
      typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
      const static int MatEledim;

      void setup(VSystem *Ptr,int _A)
      {
	A=_A;
	pSystem=Ptr;
	setupOrbitals();
      	setupTwoBodyStates();
	//	setupTwoBodyMat();
	beta=0.;
      }
      void setupOrbitals();
      void setupTwoBodyStates();
      void setupTwoBodyMat();

      void setupIndiceIn2BStates();
      vector<int> IndiceIn2BStates;

      void get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle &MatEle) const;
      void get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle &MatEle) const;
      
      DataType get2B(int a,int b,int c,int d) const
      {
	TwoBodyMatEle MatEle;
	//	get2BmatAt(a,b,c,d,MatEle);
	get2BmatOnNeed(a,b,c,d,MatEle);
	if(MatEle.empty()) return 0;
	//	return MatEle[0]+(beta-1)*MatEle[3]*pSystem->hbar_omega/A-beta*MatEle[2]*pSystem->hbar_omega/A;
	//	return MatEle[0];//-MatEle[3]*pSystem->hbar_omega/A;
	//return MatEle[0]-pSystem->hbar_omega/A*MatEle[2];
	return MatEle[0]-pSystem->hbar_omega/A*MatEle[3];
      }
      DataType get0B() const
      {
	return 0;
	//return -pSystem->hbar_omega*1.5*beta;
      }
      DataType get1B(int a,int b) const
      {
	DataType temp(0);
	if( Orbitals[a].isInAGroupWith(Orbitals[b]) )
	  //	  temp+=(1.-1./A)*(1-beta)*pSystem->get1B_p2(Orders[a],Orders[b])+beta*pSystem->get1B_p2_plus_r2(Orders[a],Orders[b]);
	  //	  temp+=pSystem->get1B_p2_plus_r2(Orders[a],Orders[b]);
	  temp+=(1.-1./A)*pSystem->get1B_p2(Orders[a],Orders[b]);
	return temp;
      }
      
      double jplus(int a,int b) const;///< J+ matrix element between two orbital a and b
      double jminus(int a,int b) const;///< J+ matrix element between two orbital a and b
      
      void printOrbitals() const;
      void printTwoBodyStates() const;

      int A;
      
      VSystem * pSystem;
      vector<int> Orders;//split from which orbital in pSystem->Orbitals, the main shell of each orbital

      int maxProtons;
      int maxNeutrons;

      double beta;
    private:
      int mapab(int a,int b) const//a<b,a map from a<b matrix to one dimension vector
      {
	int totalOrbitals=Orbitals.size();
	return totalOrbitals*a-((a+3)*a)/2+b-1;
      }
    };




  
  //****************************************************************

  const size_t numBits=980;//number of bits for config., must be greater than totalOrbitals
  ///many body basis, configurations
  class Configs
  {
  public:
    typedef System_Tz_SM_M SystemType;
  Configs(const SystemType *Ptr,int _vN,int _vZ):pSystem(Ptr),maxProtons(pSystem->maxProtons), \
      maxNeutrons(pSystem->maxNeutrons),Orbitals(Ptr->Orbitals),vN(_vN),vZ(_vZ)
    {
      //check if numBit is enough
      if(pSystem->Orbitals.size() > numBits)
	{
	  std::cerr<<"numBits("<<numBits<<") "<< "must not be less than totalOrbitals("<<pSystem->Orbitals.size()<<"), please set numBits larger in sm_system.h!!!"<<std::endl;
	  exit(0);
	}
      //check if A,Z is properly set.
      if(_vZ > maxProtons || _vN > maxNeutrons)
	{
	  std::cerr<<"number of valence protons or neutrons exceed number of orbitals(maxProtons:"<<maxProtons<<",maxNeutrons:"<<maxNeutrons<<")!!!\n";
	  exit(0);
	}
    }

    void setup(int MM,int Par,DataType Emax);///< setup configs which has total angular momentum projection MM/2 and parity Par. Par=0 for +, 1 for -
    
    vector< bitset<numBits> > configs;///< different configurations; A basis contain a sequence of bits, 0-maxProtons-1 for protons, this is in accordance with the orders of orbitals in pSystem.

    DataType getME(int i,int j) const;///return interaction matrix element between config i and j.
    double getJ2(int i,int j) const;//return J^2 matrix element between config. i and j.
    
    const SystemType * pSystem;
    const vector<SystemType::OrbitalType> & Orbitals;
    int maxNeutrons,maxProtons;

    int vN,vZ;///< number of valence nuclei.
    int dim;
  private:
    int getMM(const bitset<numBits> & config) const;///< cal the total angular momentum projection of a config.
    int getPar(const bitset<numBits> & config) const;///< cal the parity of a config.
    DataType getE(const bitset<numBits> & config) const;///< cal the energy of a config, ie. the sum of occupied sp. energy

    DataType get0B() const
    {
      return pSystem->get0B();
    }
    DataType get1B(int a,int b) const
    {
      return pSystem->get1B(a,b);
    }
    DataType get2B(int a,int b,int c,int d) const
    {
      return pSystem->get2B(a,b,c,d);
    }
  };

  //****************************************************************



}
#endif
